# TriggerSwitchD, executes commands to activate/deactive remote switches.
# (c) Copyright 2010 Pontus Strömdahl, AdhocSkill. 
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>.


# External referenses used:
# http://www.stjarnhimlen.se (accessed 2010-03-02)
# http://users.electromagnetic.net/bu/astro/sunrise-set.php (accessed 2010-03-02)
# http://www.astro.uu.nl/~strous/AA/en/reken/zonpositie.html (accessed 2010-03-02)

require "time"

module Daylight #:nodoc:all
  def self.extended(receiver)
    days = receiver.jd
    raise "method jd on #{receiver} does not return Fixnum" unless days.kind_of? Fixnum
  end

  def sunrise(at_location)
    time_for(:sunrise,at_location)
  end

  def sunset(at_location)
    time_for(:sunset,at_location)
  end
  
  private
  
  JULIAN_DAYS_UNTIL_YEAR_2000 = 2451545.0
  UNIT_CIRCLE = 360.0
  C1 = 0.0009
  C8 = 0.0053
  C9 = 0.0069
  FOR = {:sunrise => -1, :sunset => 1}

  # Calculates declination of earth
  class DeclinationOfEarth #:nodoc:all
    Obliquity_of_the_Ecliptic_for_Earth_in_rad = 23.45 * Math::PI/180.0
    def initialize(angle_as_rad)
      @declination = Math.asin(Math.sin(angle_as_rad) * Math.sin(Obliquity_of_the_Ecliptic_for_Earth_in_rad))
    end
    def to_f
      @declination
    end
  end
  
  # Holds angle ratio per degree
  class AngleRatio #:nodoc:all
    def initialize(angle)
      @ratio = angle.to_f/360.0
    end
    def to_f
      @ratio
    end
  end
  
  #Equation of Center for Earth
  class EquationOfCenterForEarth #:nodoc:all
    COEFFICIENTS = [1.9148, 0.0200, 0.0003]
    def initialize(mean_anomaly_as_rad)
      sum = 0
      COEFFICIENTS.each_with_index do |coefficient, multiplier|
        sum += coefficient * Math.sin((multiplier + 1) * mean_anomaly_as_rad)
      end
      @center = sum
    end
    def to_f
      @center
    end
  end
  
  #Hour Angle for observer in degrees
  class HourAngleForObserver #:nodoc:all
    Solar_Disk_Diameter_and_Refraction_for_Earth_as_rad = -0.83 * Math::PI/180.0
    def initialize(north_as_rad, transit)
      declination = DeclinationOfEarth.new(transit.lamda.to_rad)
      hour_angle_numerator = Math.sin(Solar_Disk_Diameter_and_Refraction_for_Earth_as_rad) - Math.sin(north_as_rad)*Math.sin(declination)
      hour_angle_denominator = Math.cos(north_as_rad) * Math.cos(declination)
      @hour_angle = Math.acos(hour_angle_numerator/hour_angle_denominator).to_degrees
    end
    def to_f
      @hour_angle
    end
  end
  
  #Calculates new transit, by calling step a number of times to improve the approximation
  class Transit #:nodoc:all
    
    attr_reader :lamda
    def initialize(west,north,date)
      @north = north
      # @what_time = what_time
      @julian_date = JulianDateApproximation.new(west,date)
      @approximation = transit_approximation(@julian_date.to_f)
      @mean_anomaly, @lamda, @julian_date_transit_approximation_as_float = @approximation
      run
    end
    
    def to_utc(what_time)
      @julian_date.hour_angle = FOR[what_time]*Float(HourAngleForObserver.new(@north.to_rad, self))
      raw = Float(@julian_date) + (C8 * Math.sin(@mean_anomaly)) - (C9 * Math.sin(2*@lamda.to_rad))
      dt = Date.jd(raw+0.5)
      Time.httpdate(dt.strftime("%a, %d %b %Y %H:%M:%S GMT"))
    end
    
    private
    def run
      number_of_approximations_to_go = 5
      while (this_is_not_the_best_approximation?(new_approximation = transit_approximation(@julian_date_transit_approximation_as_float - JULIAN_DAYS_UNTIL_YEAR_2000)) || number_of_approximations_to_go == 0)
        @approximation = new_approximation
        @mean_anomaly, @lamda, @julian_date_transit_approximation_as_float = @approximation
        number_of_approximations_to_go -= 1
      end
    end
    
    def this_is_not_the_best_approximation?(new_approximation)
      @approximation != new_approximation
    end
    
    def transit_approximation(julian_date_transit_approximation_as_float)
      mean_anomaly = (357.5291 + 0.98560028 * julian_date_transit_approximation_as_float) % UNIT_CIRCLE
      __new_transit__(__lamda__(mean_anomaly.to_rad))
    end
    
    def __new_transit__(m_lamda)
      mean_anomaly, lamda = m_lamda
      transit = (@julian_date.to_f + (C8 * Math.sin(mean_anomaly)) - (C9 * Math.sin(2 * lamda.to_rad)))
      [mean_anomaly, lamda, transit]
    end

    def __lamda__(mean_anomaly)
      center = Float(EquationOfCenterForEarth.new(mean_anomaly))
      lamda = (mean_anomaly.to_degrees + center + 360.0/2 + 102.9372)
      [mean_anomaly, lamda % UNIT_CIRCLE]
    end
  end
    
  #Julian date approximation by observer westerly location
  class JulianDateApproximation #:nodoc:all
    attr_accessor :hour_angle
    def initialize(west,date)
      @date = date
      @west = west
      @approximation = approximate(@west,@date)
    end
    def hour_angle=(angle)
      @approximation = approximate(@west, @date, angle)
    end
    
    def to_f
      @approximation
    end
    
    private
    
    def approximate(west, date, hour_angle=0)
      number = ((date.jd - JULIAN_DAYS_UNTIL_YEAR_2000 - C1) - Float(AngleRatio.new(west))).round
      JULIAN_DAYS_UNTIL_YEAR_2000 + C1 + Float(AngleRatio.new(west + hour_angle)) + number
    end
    
  end
    
  def time_for(what_time,at_location)
    __prepare_float__
    Transit.new(at_location[:west],at_location[:north],self).to_utc(what_time)
  end

  def __prepare_float__
    Float.class_eval {include DayFloat}
    def __prepare_float__
    end
  end
  
  module DayFloat
    D180_PI = 180.0/Math::PI
    PI_D180 = Math::PI/180.0
    def to_rad
      self * PI_D180
    end
    def to_degrees
      self * D180_PI
    end
  end
  
end